The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/spline

Text File | 1995-03-31 | 13KB | 490 lines

Article 4947 of comp.sys.handhelds:
Path: en.ecn.purdue.edu!noose.ecn.purdue.edu!samsung!usc!sdd.hp.com!elroy.jpl.nasa.gov!ames!haven!umbc3!math13.math.umbc.edu!rouben
From: [e-mail] (Rouben Rostamian)
Newsgroups: comp.sys.handhelds
Subject: SPLINE v2.0 for HP48 (repost)
Message-ID: <[e-mail]>
Date: 11 Mar 91 15:16:33 GMT
Sender: [e-mail]
Reply-To: [e-mail]
Organization: Mathematics Department University of Maryland, Baltimore County
Lines: 475

I had trouble with posting this porgram earlier today.  This is a repost.
I apologize if you received duplicate copies.
--
Rouben Rostamian                          Telephone: [phone]
Department of Mathematics and Statistics  e-mail:
University of Maryland Baltimore County   bitnet: [e-mail]
Baltimore, MD 21228,  U.S.A.              internet: [e-mail]

---------------------------------------------------------------------------
This is version 2.0 of SPLINE.  SPLINE generates a piecewise cubic and twice
continuously differentiable interpolation y(x) of a set of points (x_i,y_i),
i=1,2,...,n.  It is assumed throughout that x_1 < x_2 < ... < x_n.

SPLINE's default action is to generate a _natural_ cubic spline, i.e., the
second derivative y'' vanishes at the end points x_1 and x_n.  The default
action of SPLINE may be modified by specifying optional switches which are
described later.

------- INPUT --------------------------------------------------------
SPLINE reads its input from the stack.  The n coordinate points may be
specified in two different formats:

THE COORDINATE-PAIRS INPUT FORMAT:
n+1: (x_1,y_1)
n:   (x_2,y_2)
     ...
2:   (x_n,y_n)
1:           n

THE ARRAY INPUT FORMAT:
2:   [ y_1 y_2 ... y_n ]
1:           [ x_1 x_n ]

In the ARRAY input format the interval [x_1 x_n] is automatically divided
into n equally spaced nodes.  The COORDINATE-PAIRS format is useful if
the nodes are not equally spaced.

----------- OUTPUT ----------------------------------------------------
The output of SPLINE is a *program* which can be used as a user-defined
function.  The program, which is placed on level 1 of the stack, has the
following general format:

1: << -> X
      <<
          Description of the spline curve here
      >>
   >>

This program may be stored in a variable, say TRY, and may be
evaluated as "TRY(X)" (algebraic mode) or as "X TRY" (RPN mode.)
TRY(X) may be plotted with the usual plotting commands.

-------- OPTIONAL SWITCHES ---------------------------------------------
o  SPLINE by default imposes the natural boundary conditions
y''(x_1) = y''(x_n) = 0.  It is possible to specify instead the first
derivatives a := y'(x_1) and b := y'(x_n) as boundary conditions.
For this, enter the data in one of the two formats described
before, then push the list { a b } into the level 1 of stack.

o  SPLINE can also return the first derivative y'(x) and the second
derivative y''(x) of the cubic spline interpolant.  To get these, enter
data as before, optionally enter the list { a b } from the
previous paragraph, then push the characters 'D1' into level 1
to compute y'(x), or 'D2' to compute y''(x).

---------- EXAMPLES --------------------------------------------------
Example 1:
    5:  (0,0)
    4:  (1,1)
    3:  (2,4)
    2: (4,16)
    1:     4

SPLINE returns the program:

1: << -> X
  <<
    CASE 'X' 2 >=
      THEN '2.60869565217*(4-X)^3/2/6+6.86956521739*(X-2)+2.26086956522'
      END 'X' 1 >=
      THEN '(2.34782608696*(2-X)^3+2.60869565217*(X-1)^3)/6
                     +2.95652173913*(X-1)+.608695652173'
      END  '2.34782608696X^3/6+.608695652173*X'
    END EVAL
  >>
>>

Example 2a:
    2: [ 0 1 0 1 0 1 0 ]
    1:           [ 0 9 ]

The second derivative is zero at the end points.

Example 2b:
    3: [ 0 1 0 1 0 1 0 ]
    2:           [ 0 9 ]
    1:           { 0 0 }

The first derivative is  zero at the ends.   Replacing { 0 0 } by { 1 -1 }
makes the first derivative equal 1 and -1 at the ends.

Example 2c:
    4: [ 0 1 0 1 0 1 0 ]
    3:           [ 0 9 ]
    2:           { 0 0 }
    1:              'D1'

Computes the first derivative y'(x) of the spline y(x) of example 2b.
Replacing 'D1' with 'D2' will computes the second derivative.   It is
instructive to save y(x), y'(x), and y''(x) into variables Y0, Y1, and Y2,
and PLOT { 'Y2(X)' 'Y1(X)' Y0(X)' } with XRANGE set to 0,9 and AUTO.

-------- REFERENCE -----------------------------------------------------
Stoer and Bulirsch, Numerical Analysis

-------- REMARKS -------------------------------------------------------
SPLINE does not use, change, create or modify any global variables.
It does not modify parts of the stack it does not own and does not
alter any system flags, although the calculator has to be in the
symbolic mode for SPLINE to operate.   It clears _user_ flags 6,7,8,9.

-------- NOTES ---------------------------------------------------------
SPLINE V2.0 is completely different from an earlier version (no version
number) that I posted to comp.sys.handhelds a few weeks ago.  This new
version generates code which executes about 4 times faster than the
previous version.  It also has many additional features.  I will not
describe the differences here because because I consider the previous
version obsolete.

-------- PROGRAM OBJECT CHECKSUMS --------------------------------------
Checksum: #3B69h
Bytes:      2212

-------- COMMENTED PROGRAM (Uncommented program follows) --------------

%%HP: T(3)A(D)F(.);
\<<
@ Display version while working
"     SPLINE V2.0
  " 1 DISP

  6 CF 7 CF 8 CF 9 CF        @ Prepare flags 6-9

  IF DUP TYPE 6 SAME        @ Check if 'D1' or 'D2' are specified
  THEN
    CASE DUP 'D1' SAME
      THEN 6 SF            @ Set flag 6 if 'D1' specified
      END DUP 'D2' SAME
      THEN 7 SF            @ Set flag 7 if 'D2' specified
      END DUP \->STR ": Unknown flag" + DOERR
    END
    DROP            @ Drop 'D1' or 'D2' from stack
  END

  IF DUP TYPE 5 \=/        @ See if the end-derivatives are specified
  THEN { 0 0 }            @ If not, insert {0 0} (as a place holder)
  ELSE 8 SF            @ If yes, then set flag 8
  END

  IF OVER TYPE 0 SAME        @ If we have a number n in level 2 then we
  THEN OVER 2 + ROLLD        @ expect n pairs of coordinates above it
  ELSE SWAP ROT            @ Otherwise we expect two arrays in levels 2
  END                @ and 3.  In either case, we move the { s1 s2 }
                @ list form level 1 to the top of the stack.

  IF DUP TYPE 0 \=/        @ If we don't have a number in level 1
  THEN 9 SF            @ Then the coordinates are given as arrays
  DUP SIZE EVAL            @ Determine the size of array
  END

                @ Set up the local variables
     DUP 1 - { } { } 0  0   0  0 0 0 0 0   0    0   0   0
  \->   n  k  x   y  h \Gl \Gm s d m a b \Gl0 \Gmn s11 s1n

  @ Begin the main program
  \<<

    IF 9 FC?C
    THEN 1 n            @ Flag 9 is clear so data is in coord. pairs
      FOR j            @ Convert coordinate pairs into lists x and y
          C\->R 'y' STO+ 'x' STO+
      NEXT
    ELSE
      OBJ\-> EVAL \->LIST 'y' STO    @ Store array of y_i into the list y
      OBJ\-> DROP OVER - k /        @ Compute the mesh size
      0 k
      FOR j
        DUP j * 3 PICK + 3 ROLLD    @ Generate the x mesh
      NEXT
      DROP2 n \->LIST 'x' STO        @ Store the x values into the list x
    END

    EVAL 's1n' STO 's11' STO        @ Read the end-derivative values

    x EVAL                 @ Compute the list of
    1 k                    @ interval lengths h_j = x_{j+1} - x_j
    FOR j OVER - n ROLLD
    NEXT
    DROP
    k \->LIST 'h' STO

    1 k                    @ Compute the list of slopes s_j
    FOR j                @ s_j = ( y_{j+1} - y_j ) / h_j
      y j GETI 3 ROLLD GET - NEG h j GET /
    NEXT
    k \->LIST 's' STO

    @ Compute the elements d_j of the list d:
    IF 8 FS?
    THEN                @ End-derivatives are specified
      1 '\Gl0' STO
      1 '\Gmn' STO
      s 1 GET s11 - h 1 GET / 6 *
    ELSE 0
    END

    @ Still computing d:
    1 n 2 -
    FOR j
      s j GETI 3 ROLLD GET - NEG h j GETI 3 ROLLD GET + / 6 *
    NEXT

    @ End of computation of d:
    IF 8 FS?C
    THEN
      s1n s k GET - h k GET / 6 *
    ELSE 0
    END
    n \->LIST 'd' STO


    @ Compute lambda_j:
    h OBJ\-> 1 - 1 SWAP
    FOR j
      DUP 3 PICK + / k ROLLD
    NEXT
    DROP
    n 2 - \->LIST '\Gl' STO

    @ Compute gamma_j:
    \Gl OBJ\-> 1 SWAP
    FOR j
      NEG 1 +
      n 2 - ROLL
    NEXT
    n 2 - \->LIST '\Gm' STO

    @ Compute the moments m_j:
    n IDN 2 *
    2 k
    FOR j
      j DUP 1 - 2 \->LIST
      \Gm j 1 - GET
      PUT
      j DUP 1 + 2 \->LIST
      \Gl j 1 - GET
      PUT
    NEXT
    2 \Gl0 PUT
    n SQ 1 - \Gmn PUT
    INV
    d OBJ\-> \->ARRY *
    'm' STO

    @ Compute a_j:
    1 k
    FOR j
      m j GETI 3 ROLLD GET - h j GET * 6 / s j GET +
    NEXT
    k \->LIST 'a' STO

    @ Compute b_j:
    1 k
    FOR j
      y j GET m j GET h j GET SQ * 6 / -
    NEXT
    k \->LIST 'b' STO

    @ Now we compute the individual arcs of the spline:
    CASE 6 FS?C
      THEN                @ Will compute y'(x)
        1 k
        FOR j
          m j 1 + GET 'X' x j GET - SQ *
      m j GET x j 1 + GET 'X' - SQ * -
      h j GET / 2 / a j GET +
        NEXT
      END 7 FS?C
      THEN                @ Will compute y''(x)
      1 k
      FOR j
        m j GET x j 1 + GET 'X' - *
    m j 1 + GET 'X' x j GET - * +
    h j GET / NEXT
      END
      1 k
      FOR j                @ Will compute y(x)
      m j GET x j 1 + GET 'X' - 3 ^ *
      m j 1 + GET 'X' x j GET - 3 ^ * +
      h j GET / 6 /
      a j GET 'X' x j GET - * +
      b j GET +
      NEXT
    END

    @ Create the output program:
    "\<<\-> X\<<CASE " n ROLLD
    k 2
    FOR j
      'X' " " + x j GET + " \>= THEN " + SWAP + " END " +
      j ROLLD
      -1
    STEP
   " END EVAL\>>\>>"

   @ Concatenate all parts:
   1 n
   FOR j +
   NEXT
   OBJ\->
  \>>
\>>

----------- UNCOMMENTED PROGRAM -------------------------------------------

%%HP: T(3)A(D)F(.);
\<<
"     SPLINE V2.0
  "
1 DISP 6 CF 7 CF 8
CF 9 CF
  IF DUP TYPE 6
SAME
  THEN
    CASE DUP 'D1'
SAME
      THEN 6 SF
      END DUP 'D2'
SAME
      THEN 7 SF
      END DUP \->STR
": Unknown flag" +
DOERR
    END DROP
  END
  IF DUP TYPE 5 \=/
  THEN { 0 0 }
  ELSE 8 SF
  END
  IF OVER TYPE 0
SAME
  THEN OVER 2 +
ROLLD
  ELSE SWAP ROT
  END
  IF DUP TYPE 0 \=/
  THEN 9 SF DUP
SIZE EVAL
  END DUP 1 - { } {
} 0 0 0 0 0 0 0 0 0
0 0 0 \-> n k x y h \Gl
\Gm s d m a b \Gl0 \Gmn
s11 s1n
  \<<
    IF 9 FC?C
    THEN 1 n
      FOR j C\->R 'y'
STO+ 'x' STO+
      NEXT
    ELSE OBJ\-> EVAL
\->LIST 'y' STO OBJ\->
DROP OVER - k / 0 k
      FOR j DUP j *
3 PICK + 3 ROLLD
      NEXT DROP2 n
\->LIST 'x' STO
    END EVAL 's1n'
STO 's11' STO x
EVAL 1 k
    FOR j OVER - n
ROLLD
    NEXT DROP k
\->LIST 'h' STO 1 k
    FOR j y j GETI
3 ROLLD GET - NEG h
j GET /
    NEXT k \->LIST
's' STO
    IF 8 FS?
    THEN 1 '\Gl0' STO
1 '\Gmn' STO s 1 GET
s11 - h 1 GET / 6 *
    ELSE 0
    END 1 n 2 -
    FOR j s j GETI
3 ROLLD GET - NEG h
j GETI 3 ROLLD GET
+ / 6 *
    NEXT
    IF 8 FS?C
    THEN s1n s k
GET - h k GET / 6 *
    ELSE 0
    END n \->LIST 'd'
STO h OBJ\-> 1 - 1
SWAP
    FOR j DUP 3
PICK + / k ROLLD
    NEXT DROP n 2 -
\->LIST '\Gl' STO \Gl
OBJ\-> 1 SWAP
    FOR j NEG 1 + n
2 - ROLL
    NEXT n 2 -
\->LIST '\Gm' STO n IDN
2 * 2 k
    FOR j j DUP 1 -
2 \->LIST \Gm j 1 - GET
PUT j DUP 1 + 2
\->LIST \Gl j 1 - GET
PUT
    NEXT 2 \Gl0 PUT n
SQ 1 - \Gmn PUT INV d
OBJ\-> \->ARRY * 'm'
STO 1 k
    FOR j m j GETI
3 ROLLD GET - h j
GET * 6 / s j GET +
    NEXT k \->LIST
'a' STO 1 k
    FOR j y j GET m
j GET h j GET SQ *
6 / -
    NEXT k \->LIST
'b' STO
    CASE 6 FS?C
      THEN 1 k
        FOR j m j 1
+ GET 'X' x j GET -
SQ * m j GET x j 1
+ GET 'X' - SQ * -
h j GET / 2 / a j
GET +
        NEXT
      END 7 FS?C
      THEN 1 k
        FOR j m j
GET x j 1 + GET 'X'
- * m j 1 + GET 'X'
x j GET - * + h j
GET /
        NEXT
      END 1 k
      FOR j m j GET
x j 1 + GET 'X' - 3
^ * m j 1 + GET 'X'
x j GET - 3 ^ * + h
j GET / 6 / a j GET
'X' x j GET - * + b
j GET +
      NEXT
    END
"\<<\-> X\<<CASE " n
ROLLD k 2
    FOR j 'X' " " +
x j GET +
" \>= THEN " + SWAP +
" END " + j ROLLD
-1
    STEP
" END EVAL\>>\>>" 1 n
    FOR j +
    NEXT OBJ\->
  \>>
\>>